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Abstract 

We present and analyse a space-time discontinuous Galerkin method for wave propaga¬ 
tion problems. The special feature of the scheme is that it is a Trefftz method, namely 
that trial and test functions are solution of the partial differential equation to be discret- 
ised in each element of the (space-time) mesh. The method considered is a modification 
of the discontinuous Galerkin schemes of Kretzschmar et al. [25] and of Monk and 
Richter [28]. For Maxwell’s equations in one space dimension, we prove stability of the 
method, quasi-optimality, best approximation estimates for polynomial Trefftz spaces 
and (fully explicit) error bounds with high order in the meshwidth and in the polyno¬ 
mial degree. The analysis framework also applies to scalar wave problems and Maxwell’s 
equations in higher space dimensions. Some numerical experiments demonstrate the the¬ 
oretical results proved and the faster convergence compared to the non-Trefftz version 
of the scheme. 

Keywords: Discontinuous Galerkin method, Trefftz method, space-time finite ele¬ 
ments, wave propagation, Maxwell equations, a priori error analysis, approximation 
estimates. 

AMS subject classification: 65M60, 65M15, 41A10, 41A25, 35Q61. 


1 Introduction 

In this paper we analyse the Trefftz discontinuous Galerkin (Trefftz-DG) finite element 
method for the space-time discretisation of time-dependent wave propagation problems. 

Space-time finite element approximations of time-dependent wave problems constitute 
a viable alternative to methods based on finite element discretisations in space combined 
with time-stepping schemes. They provide a setting where high-order approximation in both 
space and time is obtained by simply increasing the order of the basis functions, spectral 
convergence of the error in the space-time domain can be achieved by p- refinement, and 
the hp-refinement techniques developed for spatial discretisations can be directly imported 
to space-time meshes and polynomial spaces. In particular, when local mesh refinement 
in space-time is implemented, the smallest space elements do not put constraints on the 
CFL condition in the rest of the domain. Space-time finite elements with continuous space 
discretisations were introduced in [12,13,21,22,24,33] (several references to earlier works can 
be found in [22]); space-time DG methods have been studied in [5,10,16,26,28,34], 

In order to reduce the number of degrees of freedom needed to obtain a given accur¬ 
acy, the space-time DG approach can be combined with the use of Trefftz approximating 
spaces, namely discrete spaces whose elements are piecewise solutions of the equation to be 
discretised. 
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While the literature on Trefftz finite elements for time-harmonic wave propagation prob¬ 
lems is nowadays quite developed (see e.g. [4,6,11,14,23,29,32] for different approaches using 
Trefftz-type basis functions and e.g. [2,15,19,20] for theoretical analyses), Trefftz methods 
for time-dependent wave problems are relatively new. For time-dependent scalar wave prob¬ 
lems, a global Trefftz approach was first introduced in [27]. Space-time Trefftz methods 
enforcing the continuity constraints by Lagrange multipliers have been designed, based on 
the second order in time formulation [30,35], while a Trefftz-DG method has been introduced 
in [25] for the electromagnetic wave problem in the one-dimensional spatial case and extended 
in [8] to the multidimensional case and to transparent boundary conditions. Numerical tests 
have shown that these methods actually deliver high-order space-time convergence, but no 
theoretical analysis was available. 

Here, we focus on the discontinuous Galerkin approach, and we aim at developing an error 
analysis of space-time Trefftz-DG methods, based upon the framework developed in [19] for 
time-harmonic wave problems. 

Following [25] , the model problems we consider are the time-dependent Maxwell equations 
in one space dimension, with piecewise constant material coefficients and with either Dirichlet 
or Robin boundary conditions; they also cover the case of ID scalar wave problems formulated 
as a first order system (see Remark 2.1). We prove that the Trefftz-DG method is well-posed 
and quasi-optimal, as its bilinear form is coercive and continuous in a DG norm, that the 
L 2 norm of the error is controlled, and that the scheme is dissipative. Since the Trefftz-DG 
method we consider is not of interior penalty type, no inverse estimates are needed, thus our 
scheme and its analysis also work with non-polynomial Trefftz approximations. The analysis 
framework of section 4 immediately extends to both scalar wave problems and Maxwell’s 
equations in higher space dimensions, see Remark 3.4. What is specific to the ID spatial 
case, instead, are the best approximation properties of space-time Trefftz finite element 
spaces, and therefore the error convergence rates (see sections 5 and 6). 

For the case of polynomial approximations, we prove that the Trefftz-DG method has a 
much better asymptotic behaviour in terms of accuracy per number of degrees of freedom, for 
both the h- and the p- versions, than standard methods using full polynomial discrete spaces. 
The approximation bounds for Trefftz spaces and those for full polynomial spaces have the 
same order h p+1 , with h the meshsize and p the polynomial degree, but a local full polynomial 
space of degree p for the approximation of ( E , H) has dimension (p+ l)(p + 2) = 0(p 2 ), while 
the corresponding Trefftz space has only dimension 2p + 2 = O(p). The constants in the final 
error bounds for the Trefftz-DG method (see Theorem 6.2) are completely explicit. A further 
advantage of the use of a Trefftz method is that its variational formulation only involves 
integrals on the skeleton of the space-time mesh, thus avoiding the computational burden 
of the calculation of volume integrals. We also present some numerical tests confirming 
these theoretical results, highlighting the faster convergence (both in the meshwidth and in 
the polynomial degree) compared to the non-Trefftz version of the method, and the mild 
dependence of the numerical error on the flux stabilisation parameters, which we introduce 
for analysis purposes. 

The article is structured as follows. In section 2 we introduce the initial boundary value 
problems to be discretised and in section 3 we describe the DG formulation and its Trefftz 
version. In section 4 we prove the well-posedness of the Trefftz-DG method and the a priori 
error bounds in DG and L 2 norms. In section 5 we define local polynomial Trefftz spaces and 
prove best approximation estimates in anisotropic space-time Sobolev norms and in section 6 
we combine the previous results into hp-e rror bounds for the Trefftz-DG scheme. Finally, 
in section 7 we show the results of some numerical experiments and in section 8 we outline 
some possible extensions of the scheme and of its analysis. Appendix A contains a different 
proof of the stability bound needed to control the L 2 norm of the error; this proof gives a 
better constant than that in section 4.2 but holds only for constant coefficients. 

2 Model problems 

In this section, we introduce the model problems we are going to consider in the rest of this 
paper, namely the time-dependent Maxwell equations in one space dimension, with either 
perfectly electrically conducting or absorbing boundary conditions. We refer to Remark 3.4 
below for the three-dimensional case, while the case of the acoustic wave problem is discussed 
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in Remark 2.1 (in one space dimension) and again in Remark 3.4 (in any space dimension). 

Given a space domain 12 = (xl,xr) and a time domain I = (0, T), we set Q := fl x I. 
We denote by iiq = (rig, rig) the outward pointing unit normal vector on dQ. 

We assume the electric permittivity e = e[x) and the magnetic permeability /z = n(x) to 
satisfy e(x) > e* > 0, fi(x) > /Z* > 0, and to be piecewise constant in 12. The speed of light 
in the material is c(x) := {e{x)fi(x))~ 1 ^ 2 . 

We consider the Maxwell equations for an electromagnetic wave propagating along the 
direction x , in terms of the one-component electric and magnetic fields E = E y and H = H z , 
respectively. 

In the case of Dirichlet boundary conditions, the corresponding initial boundary value 
problem reads as 


dE d(iiH) 
dx dt 

dH d(sE) 
dx + dt 

in Q, 

(1) 

in Q, 

E(-, 0) = Eq, H(-,0) = H 0 

on 12, 


E(x l ,-) = E l , E(xr,-) = E r 

on I, 


assuming sufficiently regular data J, E 0l H 0 , El, Er. The 
perfectly electric conductor. 

case El = Er 

= 0 models a 

We also consider the case of absorbing (Robin) boundary conditions: 


f + a <f > = 0 

in Q, 


dx dt 


dH d{eE) 

in Q, 


dx dt 

(2) 

E(-,0) = E 0 , H(-,0)=H 0 

on f2, 

£ 1/2 E(x l ,-) + h 1/2 H(x l , •) = gL 

on I, 


£ 1/2 E(xr, •) - h 1/2 H(xr , •) = g R 

on I. 



Remark 2.1. The case of the scalar wave equation in one space dimension 


d 2 U . 2 d 2 U 2t 

dt 2 dx 2 

can be traced back to either (1) or (2) by setting 




and /z = 1, 


(3) 


(where c 2 = 1/e/z as above). The Dirichlet boundary conditions become ^( Xl ,■) = El 
and ^t(xR,-) = Er. The absorbing boundary conditions become e 1 ^ 2 ^- + = gL, 


0/2 §u _ 

b dt 


.. 1/2 dU _ „ 
t L dx ~ 9R- 


In the following, we introduce the space-time Trefftz-DG method and develop its analysis 
in the case of the initial boundary value problem with Dirichlet boundary conditions (1). We 
address the case of the problem with Robin boundary conditions (2) in section 4.4 and in 
Remarks 6.4 and 6.5. 


3 Space-time DG discretisation 


Let K C Q be a Lipschitz subdomain such that e and /z are constant in K; we denote by 
n k = {'R'Ki'^k ) the outward pointing unit normal vector on dK. Assume E,H £ H 1 (K). 
Multiplying the first and the second equation of (1) by the test functions vh,ve G H x {K) 
respectively, and integrating by parts we obtain 


E d - V JL 
K \ ox 


fiH 


dv R TT d'VE 


dt 


H- 


dx 


E £^£) dld( 

dt / 


+ 


J [{Evr + HvE)n x K + (hHvr + eEve)^^ ds = JJ Jvndxdt. 


(Here and in the following we omit writing the trace operator Tr : 7L 1 (AT) —► L 2 (dK).) 
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3.1 Mesh and DG notation 

We introduce a mesh Th on Q , such that its elements are rectangles with sides parallel 
to the space and time axes, and all the discontinuities of the parameters e and /r lie on 
interelement boundaries (note that the method described in this paper can be generalised to 
allow discontinuities lying inside the elements as in [25]). The mesh may have hanging nodes. 
We denote with Th = U Xg7 - h the mesh skeleton and its subsets: 

Th° T := the union of the internal horizontal element sides (t = constant), 

Th er := the union of the internal vertical element sides (x = constant), 

T° ■= [xl,x r ] x {0}, 

Tj, := [xl,x r \ x {T}, 

T r h := {x L \ x [0. T], 

:= {x R \ x [0, T], 

We define the following broken Sobolev space: 

H\Th) := {u £ L 2 (Q), v\k e H\K) MK e Th} 
and the averages and jumps of functions <f> £ H 1 (Th)'- 

m ■■= \ ^ ) ° n ( 9 A 'i n dK 2) c ^r, 

l<t>lx ■= <t>\ Kl n x Kl + 0| K2 n x K2 on (dK x (~l dIC 2 ) C T^ el , 

Ult ~ <t>\ Kl nWj + <Ak 2 on (dKx n dK 2 ) C T% or . 

We note that is the trace of (j> from the left minus that from the right; is the trace 
from lower times minus that from higher times. On T^ or , we denote by tf>~ the trace of 
<j> £ H l (J~h) taken from the adjacent element with lower time. 

3.2 DG formulation 

We introduce a (vector) finite dimensional subspace V p (Th) C tt 1 (7/i) 2 . The space-time 
DG discretisation of the initial-boundary value problem (1) consists in finding ( Eh p , i?/ lp ) £ 
V p (7/t) such that, for all K &Th and for all (ve,Vh) £ V p (7/i), it holds 


jj K ( Ehp dx +l " Hhp dt 

,rr dv E, rP dv E 

+ Hhp dx + Ehp dt 

^ dxdt 



+ J ((EhpVH + H hp v E )n x K + ( eE hp v E + fJ.H hp VH)n t K ^ ds = JJ 

Jvr da; dt. 

(4) 

The numerical fluxes Eh p and Hh p are defined on the mesh skeleton Th as 




\ E h P 

on T% 01 , 

\Hn P 

on J r j) or , 



^hp 

on T%, 

Hhp 

on T%, 


Ehp •— < 

E 0 

lE hp } + p[H hp lx 

'«■ ?{■■ a 

on 71”, 

H 0 

{•ff/ip} + otfEhpJx 

on H, 
on Tjr-, 

(5) 


E l 

on Jjf, 

Hhp &(Ehp -£^l) 

011 T^ i 



, Er 

oh Tt 

k Hhp ”h ot(Ehp Eji ) 

on T^ , 



where a £ L 00 (J r ^ er U Tfc U T^) and /? £ L 00 (J ? ; j er ) are two positive flux parameters. The 
fluxes are chosen as upwind fluxes on horizontal mesh edges, and centred fluxes with the 
addition of a penalty jump stabilisation on vertical edges; other choices are possible. 
Summing (4) over all the elements K £ Th, we obtain the variational formulation 

seek (. E hp ,H hp ) £ V p (T h ) such that 

ctDG(Ehp,Hh p ;vE,VH) = £dg(ve,vh) V(ve, vr) £ V p (7fi), (6) 
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where 


dv H 


dv H 


dv E 


9ve ' 


a D G(E hp ,H hp -VE,v H ) ■■= - ^2 JJ {. Ehp ~^ + ^ Hhp -dt~ +Hhp ~^x' +eEhp -&f) 

K&Th K 


dx d t 


+ ( £ E h Jv E \t + pH h Jv H lt)dx + / {eE hp v E + vH hp v H ) dx 

+ [ {iE hp Jlv H }x + iH hp Jlv E }x + ulEhplxlvEjx + PlHhplxlvHjx) dt 

J *v x 

+ ( - H h pV E + aE hp v E ) d t+ ( H hp v E + aE hp v E ) d t, 


Pdg(ve,v h ) ■= // Jv H dxdt+ / (eE 0 v e + pH 0 v H )dx 
JJq Jj* 

+ / E L {y H + av E )dt+ / E r (-v h + av E )dt. 

Jr* 


(7) 


Remark 3.1. In [28], a space-time DG discretisation of linear hyperbolic systems in N di¬ 
mensions is introduced. The formulation (6) with a = /3 = \ and V p (7/i) = P p (7/t) 2 , the 
space of piecewise polynomials of degree at most p on Th, is a particular case of that of [28]. 
More precisely, in the notation of [28], assume the initial boundary value problem to be 
posed in one space dimension (iV = 1), and to have a linear hyperbolic system of dimension 
m = 2 with time derivative coefficient matrix A = (J5), and space derivative coefficient 
matrix A\ = (5 J). The solution will be the vector field u = (E,H). The matrix defining 
the Dirichlet boundary conditions takes values N | jtl = ( 2 and N|jtb = (q ), and the 
space-time mesh is taken as a Cartesian mesh aligned with the space and time axes. The 

( n* n x \ 

* f : on 

n K n K ) 

the part of dK with n R = 1, M + = f([ }), M _ = t ( G 1 J 1 !) > while on the part of dK with 
= —1, M + = 5 (_\ V), M- = Ii). With these definitions, the formulation of [28] 
coincides with (6) with 'V p (J~h) = P p (7ft) 2 and a = /3 = We remark that here we are in¬ 
terested in using different approximating spaces. Since we consider only meshes aligned with 
the axis, the semi-explicit time-stepping based on macro elements described in [28, section 3] 
is not applicable in our setting. However, the flux-splitting technique described there might 
be used to generalise the DG fluxes (5) to non-Cartesian meshes. 

Remark 3.2. Once the discrete space V p (7h) has been chosen, the variational problem (6) 
can be solved as a global algebraic linear system. Alternatively, one could partition the time 
interval [0, T] into subintervals j = 1,..., n, and solve sequentially the formulation 

in every time slab fI x [tj-i,tj\ using as initial conditions the traces of the solution in the 
previous slab. In this second case, the space-time mesh needs to be aligned with the time 
slabs. The two approaches are algebraically equivalent, but the second one is in general 
computationally more advantageous. In both cases the method is implicit in time. 


3.3 Trefftz-DG formulation 

In the following we restrict the discussion to the homogeneous initial value problem, i.e. J = 0 
in (1). We define the Trefftz space 

T , T , \,-u- 1 /" 7-\2 dv E . d(nv H )_dv H d(ev E ) _ n . ,, ^ ^ -r \ 

T{T h ) Y?e,vh) € H (Tfc) , + + mall A 

If we choose V p (7/t) C T(Th), the volume term in the bilinear form (7) vanishes and the 
formulation reduces to 

seek (. E hp ,H hp ) £V P (%) such that 

aTDG{Eh P , H hp -,VE,VH) = £tdg(ve,vh) V(ve,vh) £ V p (7h), (8) 
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where 


a T DG{Eh P ,H hp -VE,v H ) ■= / (^ E hvl v E \t + ^H h Jv H jt)dx + (eE hp v E + nH hp v H ) dx 

+ / + S^plNli + ap^iplxI^Elx + P\H hp \ x \vH\x) dt 

J Er 

+ ( - H hp v E + aE hp v E ) dt+ ( H hp v E + aE hp v E ) dt, 

Jr? Jr? 


Ptdg(v E ,vh)'-= {eE 0 ve + nH 0 vu)dx 

J K 

+ / E L (v H + av E ) dt + / E R (-v H + av E ) dt. 
Jr? Jr? 


(9) 


Remark 3.3. With the choice of V p (7/t) as in section 6 below, and with the choice a = /3 = 0, 
the formulation (8) coincides with the method introduced in [25]. 

Remark 3.4. In three space dimensions, Maxwell’s equations with zero source term read: 


V x E + 


aw 

dt 


= o, 


V x H 


0(eE) 

dt 


= 0 


in Q. 


Consider a space-time mesh whose elements are Cartesian products of Lipscliitz polyhedra 
in space and time intervals. In this case, the “horizontal faces” of the mesh skeleton are the 
polyhedra across which t varies, the “vertical faces” are those across which x varies. We 
define the jumps [v] t := (v _ — v+) on horizontal faces, [v]t := x V| + n^ 2 x V| Xj 
on vertical faces, and E® := dfl x I. In the case of Dirichlet boundary conditions x E = 

n n x g(x, t) on dfl x I , with numerical fluxes similar to those in (5) (in particular with 
Eh P = S e /ipS- - and + a[E/, p ] T on J^ er ), and with the obvious 

definition of the finite dimensional Trefftz space ~V p (Th), the Trefftz-DG formulation reads 


seek (E hp , H^ p ) G V p (7/i) such that 

a TDG^‘hpi H/ip! v iJ; )=4 D DG( V E. v ff) V(vE,Vfl-) G V p (7/t), (10) 


with 

aTOci^hp, H hpi v E , v H ) := / (eE7 • [v E ] t + yS.7 ■ [v ff ]t) dx 

•'•GT 

+ / (eE/,p • v E + nU hp ■ v H ) dx 

JK 

+ [ ( - KE/jpJ • [v ff ]x + {{H^pJ • [v E ] T + a[E ftp ] T • [v s ] T + /3[H hp j T ■ [vij] T ) dS 1 

j Rr 

+ / (H^ p + a(nf 2 x E hp )) ■ (n£ x v E ) dS, 

^tdg( v e, vh) := / (eE 0 • v E + mH 0 ■ v ff ) dx + / (nf 2 x g) ■ ( - v H + a(n^ x v E )) dS, 

where (E 0 ,Ho) is the initial datum (see also [8]). 

For the homogeneous acoustic wave equation —AU + c~ 2 ^^- = 0 in any space dimension, 
setting <7 = —VI/ and v = ^r, we have 

_ dcr 1 dv . 

Vl , + _. =0 , V-<r + --=0 me, 

In the case of Dirichlet boundary conditions with datum v = c/(x, t) on dfl x I, again with 
numerical fluxes similar to those in (5), the Trefftz-DG formulation reads 

seek (vh P , cr hp ) G ~V P (%) such that 
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OTDG( v hp, <7hp-, W, t) = £tDg(™, t ) V(w, t) G Vp{%), 


( 11 ) 


with 

a TDG( v hp, ^hp', w, t) := / (c _2 %H(+%-I T I() dx + / {c~ 2 v hp w + CT hp ■ t) dx 

+ I (f l ’/>pH T !N + ■ Hn + a[t)/,J N • [w] N + ^[ff/,p]N[r] N ) dS 1 

J ^r 

+ / (rr ■ n S2 + av hp )w dS, 

:= / (c~ 2 v 0 w + <T 0 • r) dx + / g(aw - t • no) dS, 

JK 

where cr 0 = —VC/(-,0) and uo = ^(-,0) are (given) initial data. Here, the jumps are 
defined as follows: [to] t := ( w~ — w + ) and [r] t := (r“ — r + ) on horizontal faces, [w]n := 
w \k 1 + w \k 2 u k 2 and Mn := t\ Ki ■ n x Ki + t\ K 2 ■ n\ 2 on vertical faces. 

The theoretical results proved in section 4 hold true also in these cases, mutatis mutandis, 
with very similar proofs. 


4 Analysis of the Trefftz-DG method 

In this section we prove the well-posedness of the Trefftz-DG method and its quasi-optimality 
in a mesh-dependent norm (Theorem 4.4), as well as error estimates in L 2 (Q) (Corollary 4.8). 
The analysis is carried out within the framework developed in [19] for time-harmonic wave 
problems. 


4.1 Well-posedness and quasi-optimality 

We define two seminorms on tt 1 (7)i) 2 : 


111(^,^)1112,0 


,1/2 


{veIi 


s^ve 


2 1 

L 2 (J£h° r ) 2 

T 1/2 h’H]t 

2 

1 

2 

p 1/2 v H 


a 1/2 lv E ia 


2 

£ 2 PT r ) 


/3 1/2 lv H ja 


L 2 (^“) 

2 

l 2 (.^ ur?) 
2 

L 2 (^- r ) 


cA 2 ue 




( 12 ) 


111(^,^)1112,0+ 


IIK^wdlll 


DG 


V2,,- 


£ ' V 


E 


r 1/2 M 


L 2 (E™) 


a x ! 2 vh 




2 

L 2 (^“) 




2 

L 2 (^° r ) 


a” 1/2 M 


2 

i 2 (^r r ) 


The parameters e, p enter the DG and the DG + norms only through their traces on horizontal 
edges where they are continuous. 

While I]] • IIIzjg is only a seminorm in H x {fTh) 2 , it defines a norm on T(Th). 

Lemma 4.1. The seminorm ||| ■ |||ug (and thus also ||| • |||dg+ ) is a norm on the Trefftz 
space T(7h). 

Proof. It suffices to verify that if |||(ub, vh)\\\dg = 0 then v E = vh = 0. If |||(ue, ujj)|||z,g = 
0, then the pair ( v e ,Vh ) G tt 1 (D) 2 and satisfies the initial value problem (1) with Eq = 
H 0 = 0, El = En = 0 and J = 0 (this last identity follows from ( v e ,vh ) G T(%,))■ Since 
problem (1) admits a unique solution [9, section 7.2.2c], then v E = Vh = 0. □ 

Proposition 4.2 (Coercivity). For all ( ve,vh) G H l (Th) 2 the following identity holds true: 


o.dg{ve,vh',v e ,vh) = |||(ub,i;jj)|||2,g- 


(13) 
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In particular, the bilinear form arDGi .'; •) is coercive in the space T{Th) with respect to the 
DG norm, with coercivity constant equal to 1. 

Proof. Using the elementwise integration by parts in time and space 


^ ff ~K~dxdt = [ [F] t dx + [ F dx — [ F dx 

.XX J JK Ot Jjrhor Jj,T Jjr 0 


KdT h ' 

E 

K&T h 


IK 

dF 

K dx 


dxdt = f [Fj x dt+ [ Fdt- f Fdt VF G W 1 ’ 1 ^) 
Jrr r Jr? Jr? 


(14) 


and the jump identity 


v~Mt - llv% = llvf t on Vr> G H\T h ), 


(15) 


we obtain the identity in the assertion: 


<-IDg(ve,Vh',Ve,Vh ) — — 


U + (■»?,)+ £(**»*)) dxdt 


KdT h 


(14) 


+ / ( £v e l v rlt + f-tv H \[vHjt) dx + / (ev% + p.v%) dx 
+ [ ({^INx + + <4ve\1 + P{vh\ 1) dt 

J ^r 

+ ( - vhve + otv%) dt + (vhve + av%) dt 

Jr? Jr? 

[^Nt +tivfjlv H jt - ^e[u|]t - \lA v H\t J dx 
J r? ai \ / 

+ o / ( £v e + » v h) dx + - [ {ev% + pv 2 H ) dx 

2 J n 2 J^i 

+ [ ({b-Ejflb-wlz + - \veVh\x +alv E jl + PIvhJI) dt 

./TTveAv. ✓ / 


=0 


+ / av 2 E dt + av E dt 




If? 


{12 = 15) \\\{ve,vh)\\\Ig- 

The coercivity of otugO, ■) = «dg(‘, •) in T(7/ l ) follows from Lemma 4.1. □ 

Proposition 4.3 (Continuity). The following continuity bound holds: 

I utdg(E, H;v e , v h )| < 2 |||(F, F)|| | DG +1| \(v E , v h )\\\dg V(F,iJ), (v E ,v H ) G T(7h). 
Moreover, when E e = Er = 0 , it holds that 

1 / 2 , 


/ Z Z \ L/Z 

«J?)| < \/2 ^ £ 1/2 F 0 + P 1/2 H 0 L2(Jr?) ) |||(xE,X/f)||| 




DG' 


Proof. The assertions follow from the definition of the bilinear form and of the linear func¬ 
tional in (9), the norms in (12), and the Cauchy-Schwarz inequality. □ 

Theorem 4.4 (Quasi-optimality). For any finite dimensional ~V p (Th) C T(Th), the Trefftz- 
DG formulation (8) admits a unique solution ( Eh p ,Hh p ) G ~V p {7h)- Moreover, the following 
quasi-optimality bound holds: 


. inf - {ve,v h )\\\dg+- (16) 

(ve,'uh)GV p (7/i) 


\\\(E,H)-(E hp ,H hp )\\\ DG <3 












Proof. To prove uniqueness, assume that El = E r = Eq = H 0 = 0. Proposition 4.2 implies 
Eh P = Hh P = 0. Existence follows from uniqueness. For (16), the triangle inequality gives 

II \(E,H) - (E hp ,H hp )\\\ DG < HI (E,H) - (ve,vh)\\\dg + \\\{.E h p,H hp ) - {v e ,v h )\\\dg (17) 

for all (v E ,v H ) G V p (7h). Since (E hp ,H hp ) - (v E ,v H ) G V p (77i) C T(Th), Proposition 4.2, 
consistency (which follows by construction and from the consistency of the numerical fluxes), 
and Proposition 4.3 give 

II \{E hp , H hp ) — (v E , vh)\\\ 2 dg = a TDG(E — v E , H — vh\E hp — v E , H hp — vr) 

<2\\\{E,H) - (v E ,V H )\\\DG+\\\( E hp,Hhp) - {v e ,vh)\\\dg, 


which, together with (17), implies (16). □ 

Remark 4.5. If El = Er = 0, i.e. the lateral boundary conditions are homogeneous, then 
the right-hand side functional £tdg{') is continuous in DG norm, see Proposition 4.3. This, 
together with Proposition 4.2, immediately gives a stability bound on the discrete solutions 
in terms of the data, i.e. 


\\\{Eh P ,H hp )\\\ DG < V2 ( 


: 1/2 E, 




H 1 / 2 H 0 


1/2 


L 2 E° h )- 


Otherwise, we only have stability in terms of the exact solution from (16): 


\\\(E hp ,H hp )\\\ DG <ME,H)\\\ DG+ . 


The reason why a stability bound in terms of E 0: H 0 , El, Er does not hold if El,Er ^ 0 
is that, in this case, the integrals on JFfc and Jqf in the definition of £tdg{ •) (see (9)) do 
not vanish and bounding them by the Cauchy-Schwarz inequality generates terms with vr 
on U Jjf, which only allow to bound \£tdg(v e ,vr)\ from above with |||(iie,?Ih-)|||.dg+> 
instead of the weaker norm |||(r’.E,i’ir)|||.DG : 


\£tdg(v e ,vr )| < a/2 ^ 


( 

e 1/2 F 0 

2 

+ 

v 1/2 h 0 

V 


L 2 E°) 



1/2 E r 




l 2 E°) 
1/2 E r " 2 




1/2 


\\\{ve,v h )\\\dg+- 


Remark 4.6. Let us fix 0 < £ < 1. In the definition (5) of the numerical fluxes on •7q( er we may 
substitute to the averages ^Eh p J and § Hh p J the weighted averages ^Eh p J^ and 
respectively, where we have set § := El + (1 — Q(f>\K 2 on (dKi D dK 2 ) C Ef er ■ All the 
results obtained in this and the following sections remain valid in this case. 


4.2 Estimates in L 2 (Q) norm 

By virtue of Theorem 4.4, we can control the Trefftz-DG error in DG norm; it is of course 
desirable to prove a bound on the error measured in a mesh-independent norm. Following 
the argument developed for time-harmonic problems in [29, Theorem 3.1] (see also [2, The¬ 
orem 4.1], [19, Lemma 3.7]), in Proposition 4.7 we prove that the L 2 (Q) norm of any Trefftz 
function is bounded by its DG norm, thus the error estimate in L 2 (Q) norm readily follows, 
see Corollary 4.8. 

The application of the technique of [29, Theorem 3.1] relies on a certain stability estimate 
for the following auxiliary problem: 


dv E 

d(nv H ) 


in Q , 

dx 

dt 

= <\> 

dv H 
dx + 

d(ev E ) 
dt 

= Ip 

in Q, 

v E (-, 0) 

= 0, : 

vh(-, 0) = 0 

on f2, 

ve{xl, 

-) = 0, 

v E {x R ,-) = 0 

on /, 


( 18 ) 
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for <f>, if £ L 2 (Q). More precisely, we will need a bound on the L 2 norm of the traces of Ve 
and vh on horizontal and vertical segments in terms of the L 2 (Q) norm of 


S^VE 

2 

+ 

t i 1/2 v h 





L2(F£ o.'UJPT) 


+ 

P~ 1/2 ve 

2 

+ 

a x ! 2 vh 



L 2 (E™) 



L 2 (^ er 


< C 2 

— '-'stab 


d/2, 


L 2 (Q) 




LHQ) 


V(0,V-)ei 2 (Q) 2 , (19) 


for some C s t a b > 0. We have inserted the numerical flux parameters within the third and 
fourth term on the left-hand side of (19) because this is what we need in the proof of 
Proposition 4.7 below; then the constant C s tab will also depend on a and /?. 


Proposition 4.7. Assume that the estimate (19) holds true for ( ve,vh) solution of prob¬ 
lem (18). Then, for any Trefftz function ( we,wh) £ T(7 h), the L 2 (Q) norm is bounded by 
the DG norm: 


£t 1/2 w e 


l 2 (Q) 


e x ! 2 wn 


) < Cstab II |(u>E,tUff)||| DG, 


L 2 (Q) 


with C s tab as in (19). 

Proof. Let {ve, vh) be the solution of the auxiliary problem (18). The space-time vector field 
{ve,HVh) belongs to H(div x j; Q), thus it has vanishing normal jumps across any smooth 
curve lying in the interior of Q ; in particular = 0. Similarly, ( vh,£Ve ) £ 

H{div Xtt ', Q) implies = [efelt = 0. Multiplying the functions to be bounded with the 

source terms of problem (18) and integrating by parts over the elements, we have 


JJ (WElf + WH<t >) dxdt 


= £ 

KdT h 


= -£ 
K&T h 


K 


dv H , d{ev E ) , dv E , d(iiv H )\ 

WE SI +WB ~er +w »-^ +w » ST) 


dx dt 


(d{ew E ) , dw H , dw E , d{iiw H ) , , 

-ve d —^— xe + — — v h H- tt, — vr dxdt 


K 


V dt 


dx 


dx 


dt 


=o 


=0 


+ / (ew E VEnK + tlWHVH^x + w E VHn x K + w H VEn x K ) ds 
JdK ) 

/ \ew E VE + fVWHVHjt dx 


=slwE}tVE+nl w HjtVH 


+ / (eweve + P-vjhvh) dx — ( ewe ve + tiwn vh ) dx 

+ / {weVh + w H v E jx dt 


=0 


=0 


= [we] x vh + I»I ] VE 


- / ( w E VH + w h Ve )dt+ / ( w E VH + w h v e ) dt 

JTf J Ef 

h =0 =0 

— 111 {vJe , VJh ) 111 DG 

•(2 / (ex| + pv 2 H ) dx + / (/3 _1 x| + a~ x v 2 H ) dt + / a _1 x| f dt) 

V JEh or UEh JE7 e r J Eh GEh / 


1/2 


(19) 

< \Pl C s tab 111 {vJE, VJh ) 111 DG ( h 1/2 * 
Since 


l 2 (Q) 


+ 


i x / 2 if 


L 2 (Q) 


1/2 


n 1/2 w e 


L 2 (Q) 


+ 


£ V 2 W H 


1/2 


= sup 


ff 0 (w E if + w H (f) dx dt 

V2 ’ 


L 2 (Q) ) (0,,/<)££ 2 (Q) 2 (||e 1/2 <( ) || i 2 ( Q ) + 11 d 1/2 V> 11^2(0)) 


\L 2 (Q). 
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we obtain the desired estimate. 


□ 


Recalling that the error ((E — Eh p ), (H — Hh P )) G T(Th), and combining Proposition 4.7 
and the quasi-optimality in DG norm proved in Theorem 4.4, we obtain the following bound 
on the Trefftz-DG error measured in L 2 (Q) norm. 

Corollary 4.8 (Quasi-optimality in L 2 (Q)). Under the assumptions of Proposition f.l, for 
any finite dimensional Trefftz space V p (T h) C T(7h), the solution (Eh p ,Hh p ) of the Trefftz- 
DG formulation (8) satisfies the bound 


( /j~ 1 ^ 2 (E — E hp ) 2 + e-V 2 {H-H hp ) 2 ) 

V L 2 (Q) L 2 (Q) / 


1/2 


L2 (Q) ■ 

fs 3 a/ 2 C sta b inf \\\(E,H) - (ve,v h )\\\ DG+i 

(ve ,Vh)&V p(Th) 


with C s tab as in (19). 

We prove now the stability bound (19) for the solution of the initial auxiliary problem 
(18), with additional assumptions on the flux parameters a and /3. The proof is based on 
differentiation of an energy functional and Gronwall’s lemma. In case of constant material 
parameters e and /i, one can derive the bound (19), based on an exact representation of the 
solution of (18), with a better constant C s t a b than that of Lemma 4.9, with no additional 
assumptions on a and /3, see appendix A; the generalisation of this argument to higher space 
dimensions and general geometries is not straightforward. 

We introduce some notation. For an element K £ 7Ji, we denote by h X K its horizontal 
edge length and set h x := max^gr h h x K . For a face / C Jfff* , / = dK\ U dK 2 , we define 

h x f := min {h x Kl ,h x K2 }, £/ := max{£ A 'i,£A' 2 }, Uf ■= max{/.t^ 1 , Uk 2 }, 

while for a face / C Tfc U Jqf, / C dK, 


hf := h x K , £f := Ek, 


Uf ■= UK- 


Lemma 4.9. Assume that the flux parameters a and /? have the following expressions on 
any face fC Ff er U U T*: 


h x 

“1/ = a l“7 £ /< 
n f 


P\ f 



( 20 ) 


where a and b are positive constants independent of the mesh size, the material coefficients, 
and the local approximating spaces. 

The solution (v e ,vh) of the initial auxiliary problem (18) satisfies the stability bound 
(19) with 


Cstab — Wic 


T 2 

■e 


i{a-\b^} 


/ 4T^ 
\ h x 


+ TT C oo + 2c oo A ? hore 


( 21 ) 


where we have set iVhor := #{L such that (x,t) £ U Jyf for some Xl < x < and 
Coo ■= ll C ll_L=o(Q)- 

Proof. We assume that <p and if are continuous and compactly supported in Q\ the general 
case will follow by a density argument. 

Let ( ve,vh ) be the solution of problem (18), and set 

£{t) ■.= - f (ev% + fiv%)dx. (22) 

^ iflxjij 


The initial conditions give £ (0) = 0, while the equations and the boundary conditions imply 


d_ 

dt 


£{£)= f fv E - 

JQx{t} ' 


d(ev E ) , d(uv H )\ 
dt ) 


dt 


+ v H 


dx = I ( - -e-(v e vh) + v E if + v H (f] da; 

v dx ) 


= -V E {x R ,t)v H {xR,t) +VE{XL,t)v H (xL,t) + / (v E lf + V H <f) dx, 

' -v-' '-v-' Jnx{t} v 7 

=0 =0 
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which in turns implies 


£(t)=£( 0 ) + ff (ve^P + chds, 


=o 


(23) 


and 


d_ 

dt 


£(t)< \ [ (sv 2 e + e V 2 + m*4 + m V 2 ) 

1 Jnxttt v J 


/f2x{t} 
= £(t) + o 


,- 1/2 


*l> 


L 2 (f!x{(}) 


+ 




da; 

2 

L 2 (Qx{t}) 


(24) 


Gronwall’s lemma as in [9, p. 624], rf(t) < a(t)r)(t)+b(t ) => r]{t) < e$a “( s ) ds ( 77 ( 0 )+ f* b(s) ds), 
applied to (24) gives 


£(t) < e*(f( 0 ) + 


,-V 2 , 


L 2 (f!x(0,i)) 




i 2 (nx(o,t)) 


(25) 


=0 


Taking into account the definition of £(£) in (22), the bound (25) allows to control the 
terms on horizontal faces. We denote by Th or the set { t € (0, T], such that ( x,t ) € ■7 r ]) or U 
for some a;L < x < a;/?}; recall that iVhor = #Ih or - Using l/eji = c 2 , we have 




L 2 Pt“U^D 

* E e*' 

tj 6 Thor 


H 1/2 v h 


L 2 (^or u ^T ) 

M“ 1/2 ^ 2 


L 2 (f2x(0,tj)) 
2 

LHQ) 




L 2 (f2x(0,tj)) 


(26) 


2 

L 2 (Q) 


< (AW-e 7 ^)^ e 1/2 </> 

Integrating (22) in (0,f) gives 

2 r 4 

= 2 / £(s) ds 
L 2 (fix(0,t)) 7 0 

(veiI> + t/>) da; dr j ds 


e 1/2 r B 

2 

+ 



L 2 (Qx(0,t)) 



( L 3) 2 


< t 


£ 1 / 2 r £ 


nx(o,s) 
2 


ce 1 / 2 ^) 


L 2 (fix(0,i)) 

2 


H 1,2 V H 


L 2 (Qx(0,i)) 


+ 


;^ 1 / 2 V’ 


i 2 (f!x(0,t)) 

2 


y/2 


L 2 (nx(o,t)) 


1/2 


which implies 


£ 1/2 U E 

2 

+ 

H 1/2 v h 


L 2 (Qx(0,£)) 



L 2 (f2x(0,t)) 


<i 2 4( 


e 1/,2 </> 


i 2 (nx( 0 ,i)) 


M 1/2 ^ 


L 2 (Jlx(0,i)) . 

(27) 


We proceed now by bounding the terms on vertical edges. We denote by a;^ the midpoint 
of K in a; direction, so that 2(a; — Xk) < h x for all (x, t) £ K , by dI\ WE and 9/v SN the union 
of the vertical and horizontal edges of K , respectively. Taking into account the expression of 
a and /? in (20) and defining for brevity A := min{a , b -1 }, we have 


/? 1/2 v e 

h x 

<A^- 
~ h x 


h x 9 
= \ n K Z 

Ux h x .. 
n n K J JK 


2 

+ 

a 1//2 i>h 

L 2 (aic WE ) 


fi~ 1/2 V E 

2 

+ 

L 2 (9if WE ) 


L 2 (dK WE) 

e~ 1/2 v h 


I 2 (ffiWE) 


d_ 

dx 


{x — Xk)h 1 v e + (x — xk)z 1 v E jdxdt 

^ SSk + e ~ ly H +2 ^ _ Xk ' ) (^ %Ve i^ + e ~ 1 vh iHF)) dxdt 
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(is) 2A 
h 
2 A 


, y | C C ( Q 

J^JJ ( + £~ 1 v% +2 (x - %)( - — (wrVr) + n^VECp + e^VH^ 

jj ^2fj,~ 1 v% + 2e~ 1 vj 1 + ( x - Xk ) 2 (y _ V 2 + e_1 ^ 2 )^ dx dt 

2A f 

+ = r— 2 | (x - x k )veVh | dx 

n x JdK sn 

Using 2veVh < (£ur + with weight £ = ec = (/re) -1 , we have the bound 


dx dt 


( P 1/2ve 

KGTh V 

_ 4A 
~ h x 


L 2 (dK WE ) 


a 1 ^ 2 vh 


L 2 (dK WE) 


M 1/2 ve 

2 4 A 

+ YV 

£ 1/2 V H 

2 Ah x 

+ „ 


2 Ah x 

+ „ 

e- 1 ' 2 ^ 


L 2 (Q) h°° 


L 2 (Q) 2 


L\Q) 2 



+ 2A 




c(ev E + HVfl) dx 


(27),(26) / 4J 12 


< A 


—j^x~ + yd + 2 c 3 00 N h0I e T 




2 

L 2 (Q) 


L 2 (Q) 


*i V V 


L 2 (Q) 


recalling that = ||c|| Loo (q). 

This, together with (26), gives the bound (19) with constant C 2 tab as in (21). 


□ 


In case of a tensor product mesh with all elements having horizontal edges of length h x 
and vertical edges of length h 4 = h x /c , the constant Cgtab is proportional to ( h x )~ 1 / 2 . We 
stress that we cannot expect a bound like (19) with (7 s tab independent of the meshwidth: 
indeed if the mesh is refined, say, uniformly, while the term in the brackets in the right-hand 
side of (19) is not modified, the left-hand side grows (consider e.g. the simple case <j> = fi, 
ip = 0, v E = 0, v H = t). 

One could attempt to derive the stability bound (19) by controlling with (</>, ip) either 
the i? 1 (Q) or the L°°(Q) norm of ( ve,Vh ), since both these norms would then control the 
desired mesh-skeleton norm. However, this is not possible, as the solution of problem (18) 
is in general not bounded in those norms. Consider for example the simple case with source 
<p = ip = ViX{o<x-t<£~ 1 } for l e N, in Q = (—3, 3) x (0,1) and with e = n = 1. The source 
is uniformly bounded in L 2 (Q) with respect to £, namely ||</>| I l 2 (Q) ~ IMIz, 2 (Q) - 1, but the 
solution ve = vh = t'/^X{o<x-t<e~ 1 } has a jump across x = t, so it does not belong to 
H 1 (Q), and ||v.e||roo(q) = \\vh Hl°°(q) = y/l is not uniformly bounded with respect to £ £ N. 


4.3 Energy considerations 

Define the continuous and discrete energies at a given time t £ [0,T]: 

£(t) ■= \ [ (eE 2 + nH 2 ) dx, £ hp (t) ~\ [ (eEl p + fj,Hl p ) dx. 

^ JQx{t} ^ 

Consider the case where El = Ur = 0 and J = 0; we have that the energy is preserved 
for the continuous problem. 

In fact, proceeding like in the first step of the proof of Lemma 4.9, taking into account 
the equations in (1), together with El = Er = 0 and J = 0, we have = 0, which 

implies £{t) = £ (0) for any t > 0. 

We turn now to the discrete case. In order to compute £hp{T) = \ fj- T {sE 2 p + /rff 2 p ) dx, 
we consider the identity 

9 — ^TDG(Ehp , -Hfap) dTDG^Ejipi Hh pi Eh p j Hhp') — ^TDG^Ehpi U/jp) 111 (Ufap, H^p) \ | Irq, 

and obtain, by simply expanding both terms and moving to the left-hand side the term 

on^, 

£h P {T) = \ ( (eE 2 + vH 2 ) dx - i / (s(E hp - E 0 ) 2 + n(H hp - H 0 ) 2 ) dx (28) 

2 Jr o 2 Je o 
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1 

2 



(4 E hplt + viHhpft) dx 



(a [-Eft,p] ^ + /3[Eftp]^) dt 


+ f (aEhpiE^ — Eh p ) + ElH^p) dt + f (otEh p (E R — Eh p ) — E R Hhp) dt. 

J rk 


Observing the signs of the terms in this expression, we note that in the case El = E R = 0 
the method is dissipative: £h P (T) < \ /^(sEq + fiH g) dx. 


4.4 The problem with Robin boundary conditions 

So far we have studied the initial boundary value problem (1) with Dirichlet boundary con¬ 
ditions only ( E = E l / r on E^ R ). In the case of the problem (2) with Robin boundary 

conditions (e 1,/2 E — gd^ECriq = Ql/r on E^ R ) the formulation and the analysis of the 
Trefftz-DG scheme follow in the same way. We outline in this section the differences. We 
denote the modified quantities with the superscript 1Z. 

We fix a new flux parameter <5 € 7 00 (7ft U 7*) satisfying 0<5*<<S<5*<1. The 
numerical fluxes (5) on 7* and 7* are modified as 


&hp ~ 


E hp -8(E hp +^/eY/ 2 H hp -e 1/2 g L ) 
E hp - 8{E hp - 0/e) 1/2 .ffftp - £~ 1/2 gR ) 


on -Fft , 

on 7ft* 


TT1Z _ 

hp 


H hp + (1 - <5)( - (e//i) 1/2 Eftp - H hp + n 1/2 g L ) 
H hp + (1 - 8)((e/ii) 1/2 E hp - H hp - /.i _1/2 5fl) 


on 7*, 
on 7ft* 


This choice arises from imposing consistency (i.e. for Eh p = E and Hh p = H we recover 
7* = E and H n = H), and imposing that the fluxes themselves satisfy the boundary 
condition (i.e. £ 1 / 2 E'^ p — = gL/rt)- We note that now the parameter a. needs to 

be defined on 7^ er only (as opposed to on 7^ er U 7* U 7* in the case of Dirichlet boundary 
conditions). 

The bilinear forms odg(-, •) in (7) and otdg(') ■) hr (9) are modified only in the terms on 
7ft and 7*; for example, cltdg{'i') becomes 

^TDG^hpi Hhp ! Ve, Vh) 

= ■■■ + J L {~( 1 - 8)E hp v H + S(n/e) 1/2 H hp v, H - SH hp v E + (1 - 5)(e/n) 1/2 E hp v E ) dt 
+ J r ((i - 8)E hp v H + 5{n/E) 1,2 H hp v H + 8H hp v E + (1 - 5)(e/n) 1,2 E hp v E ) dt. 

Similarly, the terms on lateral sides of the linear functional Itdg(') become 


&tdg( v Ei v h) — ■ ■ ■ 


+ J ^ (fe 1/2 v h + (1 - 8)g, 1/2 v E ^gLdt 

+ J r ( - 8e~ 1/2 v H + (1 - 5)g.~ 1,2 v E S j gR dt; 


the same holds for £dg(')- We also modify the terms on 7ft U 7* in the DG norm in (12) 
as follows: 


111 (^^) 111 ^ 


... + 


(1 - 5) 1/2 (e/n) 1/4 v E 


2 


5 1/2 (/j,/e) 1/4 v H 


2 


Note that now, on 7ft and 7*, both ve and vr are controlled by the DG* norm. For this 
reason, in view of establishing a continuity property for a* £)G (-, •), we define the DG* + norm 
to be equal to the DG + norm in (12) after removing the last term on the lateral sides. 

The coercivity property of Lemma 4.1 and Proposition 4.2 holds without modifications 
for a^ DG (-, •) and ||| • IIIdg 7 ^- The continuity constant of the sesquilinear form now depends 
on the parameter <5: 


Otdg( E ’ H '’ v e,vh)\ < G* |||(E, H)\\\ dg k+ \\\(ve,v h )\\\ dg -r. 
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for all ( E,H), ( v E ,v H ) G T(%) with 

C ?:=^(l + , ”“{4r ; T^}) ■ (29) 

Note that the simplest choice S = 1/2 on U Jqf gives C^ = 2\[2. As in Theorem 4.4, the 
quasi-optinrality follows: 

\\\(E,H)-{E hp ,H hp )\\\ DG K < (1 + Cf) inf \\\(E,H) - {ve,v h )\\\dg^+- (30) 

{ve,v h )€.v p(Th) 

The linear functional fy DG (-) is now bounded in the DG TC norm, even for non-zero bound¬ 
ary conditions, thus the Trefftz-DG solution depends continuously on the problem data. This 
is a slightly stronger property than that holding in the Dirichlet case, see Remark 4.5. 

Homogeneous Robin boundary conditions, as written in (2) with g R = 9r = 0, corres¬ 
pond to absorbing materials, i.e. waves hitting the boundary are not reflected back into the 
domain Q. In the non-homogeneous case, gL and g R define the wave components entering Q 
from the sides. This is reflected in the following energy identity for the continuous problem: 


£ n {T) = £ n { 0)+ [ EH At- f EHdt 

Jr? 

= £ n (G)- f (S(^/e) 1/2 H 2 + (1 - 5)(£/h) 1/2 E 2 ) dt 

+ J ^ ( Se~ 1/2 H + (1 - 5)g.~ 1/2 E) g L dt + j 5s~ 1/2 H + (1 - 5)^~ 1/2 e) g R dt. 

The second equality is derived by splitting EH = 6EH + (1 — 5)EH, then substituting in 
the first and second term the expressions of E and H , respectively, given by the boundary 
conditions in (2). Note that the value of the right-hand side is the same for any 5 G [0,1]. 
This identity is closely replicated by the Trefftz-DG discretisation: in the evolution (28) of 
the discrete energy £ Rp , the terms on T an d jrfl are substituted by 


£%,(?) = {6{ji/e?/ 2 HZ p + (1 - 6)(e/^ 2 E 2 hp ) dt 

+ J ^(Se~ 1/2 H hp + (1 - S)g,~ 1/2 E hp S jg L dt + J ^ ( - 5e~ 1,2 H hp + (1 - 5)/j~ 1/2 E hp 'jg R dt. 

Defining a on U to be equal to (1 — d)(e//z) 1 / 2 , we note that |||(iue, wjj)|||dg < 
|||(t/;E, Wh)\\\dg ' 11 for all Trefftz functions ( we,wh) G T(Th)- This guarantees that the 
result of Proposition 4.7, namely the control of the L 2 (Q) norm with the DG norm for 
Trefftz functions, holds for the DG TC norm as well. 

Combining the results sketched in this section with the approximation bounds derived 
in section 5 (which are independent of the type of boundary conditions employed), we ob¬ 
tain convergence estimates for the Trefftz-DG scheme for problem (2); this is addressed in 
Remarks 6.4 and 6.5 below. 


5 Best approximation estimates 

5.1 Left- and right-propagating waves 

In order to approximate the solutions of the Maxwell system, we decompose them into two 
components, one propagating to the right and one to the left. This allows to represent the 
solutions in terms of two functions of one real variable. In this section we describe the relation 
between fields defined in the space-time domain and their one-dimensional representations. 
In the next section this will be used to reduce the proof of approximation estimates for Trefftz 
spaces to classical one-dimensional polynomial approximation results. 

Let D = (xq, Xi) x (to, ti) be a space-time rectangle such that e and /i are constant in it. 
In correspondence to D, we define the two intervals 

: = (^o - - ct 0 ), 

. ( O -L ) 

DJ := (x 0 + do, Xi + efi), 
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and denote their length by 


h D := xi - x 0 + c(ti - t 0 ). (32) 

Their relevance is the following: the restriction to D of the solution of a Maxwell initial value 
problem posed in K x R + will depend only on the initial conditions posed on f If, U f; see 
Figure 1. 



Figure 1: The intervals f 1 


d in (31) corresponding to the space-time rectangle D. 


Let a = ( a x ,at ) £ Nq be a multi-index; for a sufficiently smooth function v , we define 
its anisotropic derivative as 


Dfv(x,1) 


Note that, if u and w satisfy 


c“‘ 


D a v(x , t) 


1 <9l a lu(:r, f) 

dS x d^ 


u(x, t) = uq(x — ct), w(x,t) = wq(x + ct), (33) 

with uq and wq defined in flf, and respectively, then 

D“u(x,t) = (-1 {x - ct), 

D“w(x, t ) = Wq“^(x + ct). 

We define the Sobolev spaces W 3,00 (D) and H 3 C (D) as the spaces of functions whose Df 
derivatives, 0 < |a| < j, belong to L°°{D) and L 2 {D), respectively. We define the following 
seminorms: 


M W*’°°(D) '■= SU P W D >Wl°°(D) . 

I a I —3 


I Hi(D) 


E ii^ii 1 > (D) • 

l«l=i 


Note that for j = 0 they reduce to the usual L°°(D) and L 2 {D) norms and we omit the 
subscript c. On the segments the H / - , ’°°(f2^) and H 3 (12^) seminorms are defined in 
the standard way. We finally define the weighted H)(D) norm (recall the definition of hp 
in (32)) 

IMIiJJ(-D) := IMIl 2 (D) + ^D K'lffi(D) ‘ (34) 

Proposition 5.1. Assume that u(x, t) = uq(x — ct) for (x, t) £ D. Then, for j £ No, 

(i) u £ Wf’°°(D) if and only if uq € W 3 ’°°(ftf ) ), and 


(ii) if u £ W 3,oc (D), then uq £ H 3 {Q.f } ), and 

\ U o\H3(n~) — hD \u\ w j,ec^D) i 
(in) if uq £ H 3 (flf,), then u £ H 3 (D), and 

M hUd) - (j + l)min{(t 1 -t 0 ), c ,T °^ } |wolL ( n-) 
Furthermore, if j = 1, 

ll u llffi(D) - ~ ll u o|li2( n -) -+ |wo| ff i(n-) ■ 
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A similar result holds for w(x,t ) = wo(x + ct), with SlJ, instead of Q D . 
Proof. For the tF;? ,oc (D)-seminorms in (i), we have 


I wi-°°(D)= SU P \\ D ? u \\l°°(d) = SU P u^ al) (x - ct) 

M=j M=j 


L°°(D) 


= M 




For the bound of |wo| H j^ n -) in (ii), we have 


\ u o\Hpn^) ~ 

= hp \u\ w i,<x>( D } 


/ 

Uq\z) 

2 

d~ < | f Ifj | sup 

u^(z) 

2 1 1 2 
= sup sup \D“u(x,t)\ 



ze f!„ 


M =j (x,t)eD 


Consider now the H 3 C (D)-seminorm. We have 


i hud) - y 

I“l -3 


\Dfu\ z dxdt 


D 


= y jj D u o\ x ~ ct ) 

i«i=j 

=(j+i) 


Uq' 1 (x — ct) 


da; dt 


da; dt 


= (j + 1) i 


< (j + 1) min 


4l i \ 2 

u o 


to 

tl 


dt, 


, 0 ) 


L 2 (xo~ct,xi —ct) 

2 f x i i 

dt, 


f Xl 1 

„,0') 

/ - 

u 0 

>x o c 


< 


L 2 (^d) J x 0 C 

(j + 1) min {(N - to), ( Xl X °^ | |mq J) 


(a) 
o 

2 

L 2 (Q~) 


L 2 (Q" 


L 2 (x—cti ,x—do) 


da; 


dx 


from which the desired bound in (Hi) follows. Note that the two terms in the curly braces 
in the last equality are equal to each other. 

The result for w is obtained in the same way. □ 

Remark 5.2. The inequality opposite to that in item (in) of Proposition 5.1 is not true. 
For example, the functions uq^(z) = V^Xlxo-cti^o-cti+e- 1 )^) (where X(a,b) denotes the 
characteristic function of the interval (a, b)) belong to L 2 (for sufficiently large <eN, 
and ||ito,i||£ 2 (Q-) = 1- On the other hand, ug(x,t) = Uo,e(x — ct) G L 2 (D) but \\uz\\ L 2 ( D ) = 

(2c£) -1 / 2 , so no bound ||wo,t|| L 2 (Q~) < C \\ue\\ L2 ^ D ^ with C independent of ( is possible. 


5.2 Local discrete Trefftz spaces 

Given a rectangle D as above, we define the corresponding local Trefftz space as 

Any Trefftz held ( E,H) G T(Z)) can be decomposed as 


( E,H)= ^ 2 ~\/2 ’ 2 1/2 ^ ’ with u = e x I 2 E + ^l 1 / 2 H and w = — [A^ 2 H. (35) 

The waves u and w satisfy (33), i.e. they propagate in the right and the left direction, 
respectively. 

Conversely, for any u 0 G C m (fl D ) and w o G C m (fl D ) for m G No, the functions u and w 
as dehned in (33) satisfy the wave equation (3) in D, and the held ( E , H) obtained combining 
them as in (35) belongs to T(D) D C m (D) 2 . 
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This suggests a construction of discrete subspaces of T (D): given p £ No and two sets of 
p +1 linearly independent functions ={</3 q , C C m (Cl D ) and $+ ={</?,},..., <^+} C 

C m (fl^) we define the space 


'Vp(D) := span 


f/<p 0 (x-ct) p 0 {x~ct)\ /p p (x-ct ) <p p (x-ct)\ 

\\ 2ei/2 5 2/z 1 / 2 /’" ‘ ’ V 2eV2 ’ 2^/2 J’ 

/^(x + ct) ^(x + ctK /<^+(a; + ct) <^+(x + ct) 

V 2s 1 / 2 ’ 2/x 1 / 2 / ’" ‘ ’ v 2eV 2 ’ 2/x 1 / 2 


which is a subspace of T(Z?) D C m (D) 2 with dimension 2(p + 1). 

By virtue of Proposition 5.1, the approximation properties of V p (D) in T(D) only de¬ 
pend on the approximation properties of the one-dimensional functions {v^ob ■ ■ •, p±}'- for all 
(E, H) G T (D) fl W^’°°(D) 2 , defining u, w, uq and wq from ( E , H) using (35) and (33), 


inf ( e 1/2 {E-E hp ) + ^ 2 (H - H hp ) ) 

hn 6V„ ID) \ V WU (D) y W£' (D) J 


(E hp ,Hh P )eV p (D) 

(35) 

“0,p£span{v3 


(36) 


inf _ (-\u(x,i) — uo, p (x — ct)+w(x,t)-wo, p (x + ct)\ w j,o a , D , 

TVo , — t<Pp 1. 

+ i I u(x,t) - u 0tP (x - ct) - w(x,t ) -t-w 0 ,p(x-l-ct)| w j,«. (D) ^ 


i"o, P espan{v+ I ...,vj+} 


Prop. 5.1 (i) 

— _ l M o~ ' u o,p|^.°°(n“) + , , 

“0,p6span{<p 0 ,...,¥>p } wo.pGspanjipQ ,...,tp p } 


inf + | wo — Wo, 


while for all (E,H) £ T(D) D H^Df 

2 

Hi(D) 


inf ( e 1 ' 2 ^ - E hp ) * + ^[H - H hp ) 2 ) 

hJeV.(D) V F HUD) F HUD)) 


(E hp ,H hp )£V p (D) 

Prop. 5.1 (Hi) 

< 


,p\wi-°°(n+) > 


(37) 


(j + 1) min {(fi - t 0 ), —-— } 


inf_ K - u 0 , pinna-) + inf |w 0 - w 0 ,p\ Hi , n i 

■^o,p^span{(/3 0 ,...,tp p } w 0iP espa,n{ipQ ,...,ipp } 


In the following, we are going to consider polynomial bases; alternative choices, e.g. 
trigonometric functions, are possible, as suggested in [30, section 3.1]. 


5.3 Polynomial Trefftz spaces 

The most straightforward choice for the space V P (D) is to take a polynomial basis: d* - = 
>I> + = {z J JjUo- (Of course, in practical implementations of the Trefftz-DG method different 
choices of the basis for the same space might be preferred, e.g. Legendre polynomials; this 
however does not affect the approximation properties of the discrete space and the orders of 
convergence of the scheme.) In this case, the general field ( Ve,vh ) £ ~V P (D) can be written 
as 

p p 

VE(x,t) = e**/ 2 ao +e~ 1 ^ 2 ^))aj(x — cty + s~ 1 ^ 2 ^^bj(x + cty, (38) 

j =i 3 =i 

P P 

Vh{x, t ) = /tW^&o+A* -1 / 2 ^ Cbj{x — cty — pT x ^ 2 V bj(x + cty, aj,bj £ R, {x, t) £ D , 

3 =1 3 = 1 

Note that the space 'Vp(D) has dimension 2p + 2, while the full polynomial space of the same 
degree (P p ) 2 has much higher dimension (p+ l)(p + 2) but similar approximation properties 
for solutions of wave equations, as demonstrated for example in Figure 2 below. Instead of 
(P p ) 2 , tensor product polynomial spaces could be considered; the space dimension would be 
2 (p + l) 2 and the approximation rates would remain unchanged. 

We recall that our goal is to control the best approximation error in the DG + norm. 
Observing its definition (12), it is enough to control the error either in L°°(Q) or in H^(Q). 
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Following the second route, we prove simple approximation bounds in H/(Q), for a general 
rectangle D C Q, which will be chosen as a mesh element in section 6. 

We first consider approximation bounds for functions with limited Sobolev regularity; 
then, we prove approximation estimates with exponential rates in the polynomial degree for 
analytic functions. 


5.3.1 Algebraic approximation 

Classical /ip-approximation results, together with a scaling argument, give the following 
one-dimensional polynomial approximation bound (see [31, Corollary 3.15], with the norms 
defined in [31, equation (3.3.10)]): for all u £ H s+1 (a , b), s,p £ N, and s < p, 


inf 

PePP([a,i>]) 


f pjp +1) 

W-«) 2 


h 


P \\L 2 (a,b) + l U 


I m(a,b) j — 


(p + s)r 


\H“+ 1 (a,b) 


(39) 


where P p ([a, &]) denotes the space of polynomials of degree at most p in the interval [a, b]. 
Combining (39) with the bound (37) proved in the previous section, the definitions of 12)3 
(31), ho (32), u, w (35), and uq,wq (33), we have for all {E,H) £ T(7 h) D W S+1,00 (I2) 2 


(P/ip >Hhp ) 


inf ( £ 1/2 {E-E hp ) 2 + p l / 2 (H-H hp ) 2 ) 

,H hp )eV p (D)\ P H\{D) P 

4 ip — s)! u 2 s +2 ( I |2 , | 1 2 

c (p + s)! hn V |uo| ff s+1 (f2h) + ' w °\ ff s+1 (n+) ) 


< 


< 


c L + s)! h D +3 { + \ w o\ W ‘+^ 


c {p + s) 

Propjj.l (i) 4 ip — a)! 2 s+3/ 

c ip + s)\ D v 

(35) 16 (p — s)! / 




1^ 


u lw s + 1 .°°(£>) T I w Iwqf + 1 ’°°(z?) 


l / 2 E 




= (D) 


i 1/2 H 


Wc 


(40) 


D {D) 


y 


This bound can be used as best approximation estimate in the convergence analysis of the 
Trefftz-DG method; we will do this in section 6. Combining a suitable variant of (39) with 
the bounds in section 5.2, one could easily obtain similar bounds in W 3,00 (D) and H/iD)- 
however, since the inequality opposite to that of item (hi) in Proposition 5.1 does not hold 
(see Remark 5.2), we are not able to obtain H/(D) norms of (E,H) at the right-hand side 
of the approximation bounds. 


5.3.2 Exponential approximation 

The degree p of the polynomial Trefftz discrete space enters the approximation bound (40) 
through the factor (p—s)!/(p+s)!, which leads to algebraic convergence with order depending 
only on the solution regularity s. Classical polynomial approximation theory shows that, if 
E and H admit analytic extension in a complex neighbourhood of D , then exponential 
convergence in the polynomial degree p is achieved. 

Indeed, if uq and wq are analytic in the complex ellipses with foci at the extrema of 12)) 
and f2j, respectively, and sum of the semiaxes equal to pd{x i — Xq + c(t± — to))/2 for pu > 1, 
then by the classical Bernstein theorem (e.g. [7, Chapter 7, Theorem 8.1]) the exponential 
convergence rate in L°°(f2^) of the polynomial approximation of uq and wq is pd- 


inf 

Pepping) 


ll«o 


inf 11^0 — ^D,BernP]j P 

Pepp(n+) K D> 


(41) 


for some constant Cd,B ern > 0 independent of p. As in (40), combining (36) and (41), we 
obtain the following exponential approximation bound in the space-time rectangle: 


inf 

( Eh P ,Hh P )£V p (D ) 


( e^iE-Ehp) + p, l / 2 iH-H hp ) ) <C DtBein pn P . (42) 

V L°° (D) L°°(D) / 
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6 Convergence rates 

We now derive the convergence rates of the Trefftz-DG method with polynomial approxim¬ 
ating spaces 

V p (7h) = {(y E ,v H ) G L 2 {Q) 2 : (ve,vh)\ k are as in (38) with p = px }• (43) 

The two main ingredients are the quasi-optimality results investigated in section 4 and the 
best approximation bounds proved in section 5. To combine them, we need to control the 
DG + norm (12) of the approximation error with its H^(Q) norm, weighted with e and /i, to 
be able to use the bound (40). To this purpose, we define the following parameters: 


Cr- := max 


a£ 1z/» 

-{dKn(^U^UKf)) ’ 

p 

1 

"sT 

1 

L“(OR:n(X-ver u ^I. UJf -H)) 

1 

0 

°{dKn^ er ) ’ 

llr'e" 1 ! 

L~(9iCn^ er ) | 


(44) 


VAT G T h . 


For any mesh element K G Th, we denote by ft^-, hf K the lengths of its horizontal and vertical 
edges, respectively, i.e. the local meshwidths of the discretisation. 

Before stating our main convergence theorem, we prove a simple explicit trace inequality 
for functions defined on rectangles. 

Lemma 6.1. Given a space-time rectangle D = (xo,Xi) x (to,ti), denote by dD SN = 
(xq, x\) x {to, ti} and dD WE = {xo, Xi} x (to, tf) the decomposition of its boundary in opposite 
sides. For all u G H 1 (D), we have the following trace estimates: 


2 Gt 2 

IMI.L 2 (d-D SN ) — ll u lli 2 (D) ' 

2 4 2 

IMIl 2 (,9D we ) — _ II u IIl 2 (ZJ) 

xq 


1 1 — to 
2 

Xi - Xq 
2 


du 


dt 


du 


L 2 (D) 

2 


dx 


L 2 (D) 


(45) 


Proof. We assume without loss of generality that D is centred at the origin, i.e. xo = — x\ and 
to = —t\. The first bound is a simple consequence of the fundamental theorem of calculus: 


I L 2 (dDS N ) 


/ 


u 2 dx 


}x( — 



< 


L\D) 


\ L 2 (D) 


+ t\ 


du 

~dt 


L 2 (D) 


The first inequality in (45) follows from t\ = (t\ — fo)/2; the second bound is derived in a 
similar way. □ 

Theorem 6.2. For all K G Th, fix pk,sk G N with 1 < Sr- < Pk and assume that 
the restriction to K of the solution ( E,H) of the initial boundary value problem (1) with 
J = 0 belongs to W Sk+1 '°°(K). Define the discrete space V p (Th) C T(7r) as in (43), and 
let ( Ehp,Hh p ) be the solution of the corresponding Trefftz-DG variational formulation (8). 
Then, the following bound holds true: 


\\\(E,H)-(E hp ,H hp )\\\ DG 

12 ^ / h x K \ / ft* 

( 6 ( c + ^) +8Ck ( 1 + c 7 ^ 


KGTh 


1,2 le „T 


SK+? 


S K 

Pk 


-( 


1 / 2 E 


w t 


<K + 1.' 


\K) 


i 1/2 H 


w: 


sjf + l.' 


\K) 


)■ 


(46) 
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If the bound (19) holds true for the solution of the auxiliary problem (18), we also have 
the following bound in L 2 (Q): 

\ 1/2 
) / 


< 


fi- 1/2 (E-E hp ) 2 + e-V\H-H hp ) 

P L 2 (Q) 

V (e 

V K£T h V 

J c+ 7^) + 8Cif ( 1 + c 

•( 

e x ' 2 E 

w! k+1 '°°(K ) + 


2 

L 2 (Q) ■ 


1/2 




SK -\~t 


l K' / 

1/2 h 


Pk 


w; k+1 -°°(k) 


(47) 


with Cstab irom (19). 

Proof. Given an element K £ Th, we denote by dK N , dK s , 5iv w , and 5A' E its North, 
South, West and East sides, respectively, North pointing in the positive time direction, and 
set <9 AT we := dK w U dK E . 

For all ( ve,vh ) £ H 1 (Th) 2 , expanding the DG + norm (12), using the definition (44) of 
(k and the trace inequalities (45), we have the bound 

V- /I - - 2 l - - 2 

£ £ 

KdT h 


W\(vE,V H )\\\bG+ 



£ 1/2 V E 

2 1 

+ 7T 

p 1/2 v H 

\2 


L 2 (dK S) 2 



s^ve 

2 3 

+ o 

p, 1/2 v H 


L 2 (dK N ) 2 



c^^ve 


p 


- 1/2 


L 2 (dK WE ) 
2 


P 1,2 vh 


ve 


L 2 (dK s ) 

2 

L 2 (dK N ) 

2 

L 2 (8K WE \(Tf UFjf)) 
2 


( 45 ) 

£ £ 

KdT h 


—V 

/ v 


z,2 (a a- we \(j(-uj«)) 
2 


a 1|/2 Wij 


L 2 (dK WE ) 


e l/2 v E 


+ 


if 


L 2 (lf) 


p 1/2 v H 


dt 


L 2 (K) 


d(p 1/2 v H ) 


L 2 (lf) 
2 




(34) 

< 


^677 


■ CKh’k 
h 


( 

S^^ve) 

2 

d{p, 1/2 v H ) 


dx 

L 2 (K) 

dx 


^ ( 6 ( c + ^) +8C *( 1 + 4 j 

\ K K 


s^ve 


L 2 (K) 

2 

L 2 (K) 


Hl{K) 


p 1/2 vh 


2 

H\{K) 


Combining this with the quasi-optimality and the approximation results gives the error 
bound: 

(16) 

\\\{E,H)-(E hp ,H hp )\\\ DG < 3 inf \\\(E, H) - (v e ,v h )\\\ D g+ 

{Ve,vh)€-V p{Th) 


£3£ 6 (c + |) + 8 c,( 1+ c|) 

IfeTh ^ K K 


2 


( 40 ) 12 

< 


inf f £ 1 / 2 (E — ve) 

(ve ,vh) £ 

£V P (K) 

^z;(( e ( c+ k) +8( *( 1+c kd ( 


K&n 


( e 1/2 A 


^ K+1 ’”(lf) 




p 1/2 (# - Vtf) 

i/2 A^-^)!y /2 

)■ 


2 \ 1/2 
H1(K)/ 


(pat + Sic)! 


(^A‘ + 


£ \ + 2 


/ 2 ff 


lVc K + 1,a °(lf) 


The bound (46) in the assertion follows by applying to the factorial terms the upper and 
lower Stirling inequalities in the form of [1, Corollary 3.3]: for all s < p £ N 


{P~s)\ ^ 

/p — s + e 2 /2i r — 1 ^ 

v2( P -»r- 

f(l -s/p)^- 1 e 2 \ 

^ (e/2) 2s2 / P 

(p + s)! ~ 

s 

v p + s + 1/6 J 

1 (p + s)P+ s - 1 

1(1 + S /P) f+1 P 2 J 

— p2s 


<1, (s>l) 
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where the last inequality follows noting that f{x) = (1 — x) * -1 (l + x)~i~ 1 4 x e 2 ~ 2x < 1 for 
all 0 < x < 1 (which in turn can be verified by checking the convexity of log / and its limit 
values for x — > 0 and 1). 

The estimate in L 2 (Q) norm follows from combining this bound with Corollary 4.8. □ 

Remark 6.3. The constant in the first brackets at the right-hand side of the bound (46) 
controls the loss of accuracy due to anisotropic elements. If all mesh elements K satisfy 
h x = ch t , then the constants reduces to 2(3c + 4£ xY^ 2 ■ 

From bound (47) we see that, if we fix a constant polynomial degree px = p , consider 
a sequence of quasi-uniform meshes with h x ~ ch t and discretise an initial boundary value 
problem with constant coefficients and whose solution is sufficiently smooth, then the L 2 
norm of the Trefftz-DG error converges in the meshwidth with algebraic order equal to p + 1 
(note that in this case C s tab ~ h x X ^ 2 in (21)). Theorem 5.2 of [28] (see also Remark 5.2 
therein) gives the slightly lower order of convergence p + 1/2 for full polynomial spaces (thus 
with higher dimension at the same polynomial degree) in a much more general context 
(higher space dimensions, more general hyperbolic systems, non tensor-product meshes). 
The numerical experiments in [28, section 7] (confirmed by those in section 7 below) recover 
the higher order p + 1. 

The estimates of Theorem 6.2 suggest to refine the space-time mesh and reduce the 
local polynomial degree in the elements where the solution has lower regularity. The mesh 
refinement might be determined a priori knowing the locations of possible singularities in (the 
derivatives of) the initial data, discontinuities of the material parameters and non-matching 
of (the derivatives of) the initial and boundary conditions, and propagating them into Q 
along the characteristics. 

Remark 6.4. For the Robin initial boundary value problem (2) described in section 4.4, an 
analogous convergence result to Theorem 6.2 holds. After substituting the DG TC norm in 
place of the DG norm, the bounds (46) and (47) hold with a factor (1 + C ™)/3 multiplied 
to the right-hand side, with C^ as in (29) (due to the different quasi-optimality bounds in 
(16) and (30)), and with (x substituted by 


(x '■= max < II ae 11 


L°°(dK n ^ er ) 5 


| a 1 p, 1 




-i I 


L°°(dKnF% 

(i-5)M- 1/2 


I r l e~ l 


L°°(dK n(^u^«)) 


\L°°(dKDT™ T ) ’ 
L°°(dKnF™ r ) ’ 
5 { ep)~ 1 ^ 2 

L°°(dKn(^ U.F«)) 


6.1 Exponential convergence for analytic solutions 

If the solution ( E,H ) is analytic in a complex neighbourhood of some mesh elements, for 
these elements one can replace the approximation bound (40) with the exponential one (42) 
in the last step of the proof of Theorem 6.2. This suggests the design of suitable /ip-mesh 
refinements which can provide exponentially convergent Trefftz-DG discretisations. 

In the rest of this section, we discuss sufficient conditions on the problem data such that 
analyticity of the solution is guaranteed in a complex neighbourhood of all the elements in 
the mesh. 

We assume that the coefficients e and p are constant throughout Q , and, for simplicity, 
that the Dirichlet boundary conditions are homogeneous, i.e. El = Ex = 0. 

We first note that the solution (E,H) of the initial boundary value problem (1) is the 
restriction of the solution (E, H) of the similar problem posed on R x K. + with initial condi¬ 
tions l?(-,0) = Eq and H(-, 0) = H 0 , where E 0 is the 2(xr — Xi)-periodic extensions of E 0 
odd around the point Xl, and Hq is the 2(xx — a.’i)-periodic extensions of Hq even around 
the same point. 

We now assume that 

Eq and Hq are analytic in the complex strip S r := {z € C, | Imz| < r} 

for some r > 0. As in (35) and (33), we decompose the (extended) initial conditions in the 
components uq := e 1 / 2 Eq + and wq '■= e}^ 2 Eq — p}/ 2 H 0 , which are also analytic in S r . 
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(What we actually need is only that uq and wq are analytic in a sufficiently large complex 
neighbourhood of the finite segments Qq and fig, respectively.) 

For every mesh element K as above, we fix hx '■= length(fi^). The complex ellipses with 
foci at the extrema of f2^- and sum of the semiaxes equal to r + \/r 2 + h ? K /4 are contained 
in the strip S r . So the exponential approximation bound (42) holds with D chosen as K and 
exponential rate 

p K ■= Zr/hK + yj 1 + 4 r 2 /h 2 K > 1. (48) 

Proceeding as in the proof of Theorem 6.2, we see that the solution of the Trefftz-DG 
discretisation converges with exponential rates: 


\\\(E,H)-{E hp ,H hp )\\\ DG <$y/2 J2 (h x K +2( K h t K ) 1/2 C K , Bein p- 1 


KeT h 


( 2 + e~V\H-H hp ) 2 ) 

V L 2 (Q) L 2 (Q) / 


,1/2 


L 2 (Q) 

^ 6 Cstab 'y ' {h'x "p ^ Ck, Bern Px >K j 

KeTh 


(49) 


where Cx.Bem, C s t a b> Pk, (k were defined in (41), (19), (48), and (44) respectively. In par¬ 
ticular, since we have taken constant parameters, C s t a b satisfies (52) in appendix A. 

The bounds (49) ensure that, under the regularity assumptions stipulated in this section, 
the Trefftz-DG method converges exponentially in the total number of degrees of freedom on 
any fixed mesh, when the polynomial degrees px are increased uniformly. On the contrary, 
standard space-time DG methods converge as negative exponentials of the square root of 
the total number of degrees of freedom, as demonstrated e.g. in the numerical example in 
Figure 2. 

Remark 6.5. In the case of the Robin initial boundary value problem (2), convergence bounds 
similar to (49) can be proved if the functions 


v£(x) 

w£(x) 


9l({ x L - x)/c) 
£ 1 / 2 Eo(x) + p}/ 2 H 0 (x) 

£ 1/2 E 0 (x) - p} /2 H 0 {x) 
g R ((x - x R )/c) 


in (x L - cT , xl\, 
in ( x L ,x R ), 

in (x L ,x R ), 
in [x R ,x R + cT ) 


can be extended analytically in S r (or in a sufficiently large complex neighbourhood of their 
domains of definition). Note that, not only the data g R ,g Rl Eo , and Hq must be analytic, 
but their derivatives also need to match appropriately at the points (xl, 0) and (x R , 0). 


7 Numerical experiments 

In this section we present numerical results supporting the theoretical findings of the preced¬ 
ing sections. In particular, we present convergence orders for the h- and p-versions obtained 
from a series of numerical experiments. These are compared for the choice of a Trefftz basis 
and a non-Trefftz basis. For the former we employ polynomials of the Trefftz space introduced 
in (38), for the latter we choose a tensor product of Legendre polynomials in the spatial and 
temporal variables, respectively. 

7.1 Numerical method 

For all numerical experiments we employ formulation (9) with the flux stabilisation para¬ 
meters in (5) chosen asa = /3 = l/2 unless stated otherwise. In matrix form the resulting 
numerical scheme reads 

Af n+1 = Rf”, (50) 

where f" is the vector of numerical degrees of freedom at time step n, and A and R are 
assembled from terms of the bilinear form (9) acting on f at step n + 1 and n , respectively. 
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Following (38) the Trefftz basis consists of transport polynomials. As non-Trefftz basis 
functions we take 

v e’h OM) = L jx (x)L jt (t), 

where Lj[x ) denotes a Legendre polynomial of degree j. The polynomial degrees in space j x 
and time jt are independent of one another and chosen such that j x + jt < j • The dimensions 
of the Trefftz and the full polynomial space are given in section 5.3. The choice of complete 
polynomials is motivated by the observation that for transport polynomials the degrees in 
the spatial and temporal variable add up to j as well, and by the fact that for DG schemes 
complete polynomial spaces deliver similar accuracy as tensor product ones with less degrees 
of freedom, as demonstrated in [3]. 

7.2 Test problem 

In the following, we consider meshes composed by space-time squares of uniform sizes, parti¬ 
tioning the (1+1 (-dimensional domain Q = [0, 60] 2 . We enforce Dirichlet boundary conditions 
representing a perfect electrical conductor (PEC), i.e. El = Er = 0 in (1). 

As initial conditions we choose 

(h°o) = (l) exp (“ ( x ~ 10 ) 2 / 10 ) 5 

corresponding to a wave packet of Gaussian shape propagating in positive direction in free 
space. The wave is reflected once by the domain boundary. We normalise to one the permit¬ 
tivity e, the permeability p, and thus the speed of light c. 

7.3 Error convergence 

We consider the relative error computed in the L 2 (Q) norm on the whole space-time domain 
*=(//((*- E hp ) 2 + (H — H hp ) 2 ) dx d t j JJ^ (. E 2 + H 2 ) dx d?j ! . (51) 

We first consider convergence of the p- version for both Trefftz and non-Trefftz basis functions. 
The mesh step sizes are h x = ht = 1. In Figure 2 (left) we display the global relative error 
6q in dependence of the polynomial degree p. In both cases spectral convergence is observed. 
However, the Trefftz space of order p has a smaller dimension and from bound (49) we expect 
the error to decrease exponentially in the number of degrees of freedom in the Trefftz case, but 
exponentially only in the square root of the number of degrees of freedom in the non-Trefftz 
case. This is observed in the middle and right panels of Figure 2. 
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Ph 

JD 
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Figure 2: Convergence of the p- version. Comparison of the global relative L 2 (Q) error (51) 
for Trefftz and non-Trefftz basis functions for the propagation of a smooth wave packet in 
(l+l)-dimensional space against the polynomial degree (left), the square root of the number 
of degrees of freedom per element (middle), and the number of degrees of freedom per element 
(right), respectively. 
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Figure 3 shows convergence of the h-ve rsion for degree zero through three. Solid lines 
correspond to results obtained with the Trefftz basis whereas the dashed lines were ob¬ 
tained using the non-Trefftz basis. Uniform mesh step sizes are applied by reducing h x and 
ht simultaneously. The Trefftz method exhibits optimal algebraic convergence rates h p+1 . 
However, in the non-Trefftz case, the results seem to suggest an odd-even pattern of the 
convergence rates, with convergence being suboptimal for odd degrees (by one order). Nu¬ 
merical odd-even effects in the convergence rates of DG methods have also been reported, 
e.g. in [18, section 6.5], although it has been shown in [17] that in some situations this might 
be a mesh effect, the convergence being always suboptimal on particular meshes. 



Figure 3: Convergence of the h- version. Comparison of the global relative L 2 (Q) error (51) 
for Trefftz and non-Trefftz basis functions for the propagation of a smooth wave packet in 
(l-fl)-dimensional space vs. the mesh step size. The spatial and temporal step size h x and h t 
are decreased simultaneously. Optimal convergence is obtained with the Trefftz basis. In the 
non-Trefftz case, suboptimal orders of convergence are observed for odd polynomial degrees. 


7.4 Stability 

As the space-time Trefftz-DG method is implicit in time, a linear system of equations has to 
be solved for advancing the solution in every time step. In this regard, we investigated the 
conditioning of the Trefftz and non-Trefftz system matrices. Figure 4 shows that the increase 
of the conditioning with the polynomial degree in the Trefftz case is very mild, compared to 
the non-Trefftz case. 

The update matrix is obtained from (50) as U = A -1 R. Note that this matrix is usually 
not explicitly assembled but we solve for the system (50). In Figure 5 we show eigenvalues 
of the update matrices with Trefftz basis for degrees p = 0,..., 5. As expected from the 
stability analysis of section 4.3, all eigenvalues are on or within the unit circle. Figure 5 
numerically confirms stability, and it shows the method to be dissipative. 

7.5 Flux stabilisation parameters 

Let us shortly comment on the choice of the flux stabilisation parameters a, /3. The choice 
a = P = 1/2 corresponds to a full upwind flux [28] and a = /? = 0 (see [25] and also [26]) 
corresponds to a centred flux. In order to determine the influence onto the numerical error, 
we varied a, fi in steps of 0.1 in the range [0,1] and repeated the above example maintaining 
the same regular mesh with h x = h t = 1 and polynomial degree p = 2. Figure 6 depicts 
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Figure 4: Comparison of the condition number of the update matrix in the Trefftz and 
non-Trefftz case in dependence of the polynomial degree. 

the respective global relative error cq. The minimum error was obtained for a = 1/2 and 
/3 = 0. However, the overall variation of the error is within a factor of two. Similar results 
are obtained for different settings of the mesh step size and the polynomial degree. 

8 Conclusions and extensions 

In this paper we have analysed a space-time Trefftz discontinuous Galerkin method for linear 
wave propagation problems. The scheme corresponds to that proposed in [25], with the 
addition of jump penalisation terms. In one space dimension we proved that the formulation 
is well-posed and dissipative, and, for polynomial Trefftz trial spaces, we derived a priori hp- 
convergence bounds for its error in DG and L 2 norm. Numerical examples show the viability 
of the scheme and confirm the orders of converge. 

Several extensions of the analysis developed here can be envisaged. In higher space di¬ 
mensions, the well-posedness and the abstract error analysis in DG norm are straightforward. 
In order to prove orders of convergence for the scheme in higher dimensions, new best approx¬ 
imation bounds for Trefftz spaces must be developed, as those derived in section 5 rely on 
the exact representation of the PDE solution in terms of left- and right-propagating waves, 
which is a one-dimensional result. 

Other topics deserving further investigation are the generalisation of the analysis to 
unstructured meshes (see [5,28]) and the derivation of approximation estimates for non¬ 
polynomial Trefftz bases. In addition, it would be interesting to analyse the non-penalised 
case with a = /3 = 0, as well as other possible less dissipative (or non dissipative) variants. 
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A Stability bound in the case of constant coefficients 

In the special case of constant material parameters £ and /r, the stability bound (19) can 
be derived differently from Lemma 4.9, namely using an exact representation of the solution 
of (18). This results in a simpler constant C s t a b with linear, as opposed to exponential, 
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Figure 5: Eigenvalues of update matrices with Trefftz basis and flux stabilisation parameters 
a = /3 = 0.5. The grey solid line is the unit circle in the complex plane. 



Figure 6: Comparison of the global relative L 2 (Q) error (51) for the setup described in section 
7.2, under variation of the stabilisation parameters a and /3 for a Trefftz basis of degree two 
(see section 7.5). 


dependence on the total time T ; no additional assumptions on the flux parameters a and /3 
are required. 

Lemma A.l. Assume that e and p, are constant inCl (and thus in Q). The solution (ve, vh) 
of the initial auxiliary problem (18) satisfies the stability bound (19) with 

^stab — 4Tc(dV hor + rp/N vei ) , (52) 

where 

IVhor := ff{t, such that ( x,t ) £ J r ^' or UJ r )f for some Xl < x < xr}, 

-/V ve r := ff{x, such that (x, t) £ jFf l er U U Jqf for some 0 < t < T}, 
rj :=\cT/(x r - x L J], 
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7 := max { || (/3e) 1 + (a/x) 1 || ioo( ^ er) ; || (a M ) 1 }■ 

Proof. We assume that (f and if are continuous in Q ; the general case will follow by a density 
argument. 

First, we extend the initial problem to the entire space R. Define ve,vr,J>, if in R_x ®- + as 
the 2 (xr — ZzJ-periodic functions in x that satisfy ve\q = ve,vr\q = vr, 4 >\q = <f, if\ Q = if 
and such that ve and if are odd around xl (and consequently also around xr), and vr and 
<p are even around the same points, i.e. 


Ve(xl +x,t) = - v e {xl - X, t), v H {x L + X, t) = Vr(x L -X,t), 

I f(XL + x,t) = (f(XL — x,t), lf(XL + x, t) = — lf(XL ~ X, t), V(l,()gKxl + . 

(Note that the absolute values are (xr — a'cj-periodic in x .) Since time derivatives preserve 
parities and space derivatives swap them, the extended functions ve and vr are continuous 
and satisfy the extended initial problem 


dv E , d(nv H ) t 

-§^ + -§r = * 

dv H , d(ev E ) r 
v E (-,0)=0, v H {-,0) = 0 


in JK x 

in R x 
on R. 


Second, we split the right- and the left-propagating components. Define 

u + w 


u := eJ! 2 Ve + /jJ^vr, w := e^^Ve — \jJ^ 2 vr 


so that 


v E = 


2 e !/2 ’ 


Vr = 


2 pV 2 ■ 


They satisfy the inhomogeneous transport equations in R x R + 


du d(c 1 u) 1/2 ~ , 1/27 r dw 5 ( c lu ’) 1/27 1/2 7 

ai + S^ = e * + P *=■>' -di~— 5 ^ = e *='- 9 ' 

recalling that ( e/i ) 1 ^ 2 = c _1 , so they can be written explicitly with the following represent¬ 
ation formula (e.g. [9, section 2.1.2, equation (5)], recall that from the assumptions made in 
the proof, / and g are piecewise continuous) 

u(x,t) = / cf{x + c(s — t), s) ds, w(x,t) = — / cg(x — c(s — t),s) ds. 

Jo Jo 

We first bound the L 2 norm of u and w on horizontal and vertical segments with the data 
f,g ; from the triangle inequality ( ve,vh ) will be bounded by <f and if, and the bound for 
Ve and vr will follow. For all 0 < t < T 


r x R / P \ 2 

ll u ('i Oil L 2 (Q) = / (/ cf(x + c(s-t),s)dsj dx 

J XL ' J o ' 

f XR /'* 9 

<tc 2 / \f(x + c(s — t), s) |“ dsdx 

Jxr. Jo 


I XL Jo 

ft fXR+c(s — t ) 


= tp 


< 2 tc‘ 


= 2 t.c 2 


/o Jxl~\-c(s — t) 
ft pXR + c(s — t ) 


\f(y,s)\~dy ds 


If 

Jo Jxr 


0 j XL~\-c(s — t) 
2 


(e\(f{y, s)| 2 + ii\if{y, s)| 2 ) dyds 


z 1 ' 2 cf 


L 2 (nx(o,t)) 


m 1 / 2 V> 


L 2 ( f 2 x (0, t )) 


(the last equality follows from the symmetries of cf and if which ensure the equality of 
their L 2 norms on the rectangle ( xl,xr ) x (0 ,t) and on the parallelogram with vertices 
(, Xl — ct, 0), (xr — ct, 0), ( XR,t ), ( XL,t )). Similarly, for all ie!l 


t(a 


Il 2 (o,t) 


/ (/ c f( x + c (s-i),s)dsj 


d t 
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nT rt 

< c 2 / t / \f{x + c(s — £), s) | 2 ds d£ 

Jo Jo 

n T 

t\f(x + c(s — t), s) | 2 dt ds 


= c 


= c 


/ / (-—- +s)|/( 2 /,s)| 2 -dyds 

Jo J x — c (T—s) v c ' c 


'0 J x—c(T—s) 

< 2Tc /" [ (e\<j)(y, s)| 2 + /i|^(y, s)| 2 ) dy ds 

./O Jx—c(T—s) 


< 2 Tc 

< 2Tc 


£ 1 / 2 (j) 
cT 

XR - X L 


L 2 ((x—cT,x) X (0,T)) 
2 


r 1 / 2 ^ 


L 2 (Q) 




L 2 ((x—cT,x) X (0,T)) 


L 2 (Q) 


where [cT/(a;fi — £&)] is the number of times Q must be replicated in the x direction to cover 
the left domain of dependence of {x} x (0, T). 

Analogous bounds can be proved for the left-propagating term w. 

Going back to the solution of the auxiliary problem, we obtain (using ev\ < (u 2 + w 2 )/ 2, 
fivjj < (u 2 + w 2 )/ 2 and summing over all vertical and horizontal segments) 


£ 1/2 v e 


° r U^) 


+ 




+ 

0 1/2 V E 

2 

+ 

P 

1 

^0 

c* 

UJ 



LHJTT) 

1 II 


L 2 (^”U^) 

2 

L 2 (^ er U/(UJ( l ) 


7 1 / 2 m 


h h h - 

2 


— + H U; llL 2 (^ or U^) + 

< 4Tc(cN hor + \cT/(x R - ZzJllAver) ^ e 1/2 J> 


L 2 (^T r U^U^«) 

+ ^ 

l 2 (Q) 


7 1 / 2 W 


2 

L 2 (Q) 


L 2 (^ve r u^u^«) 
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